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Abstract 

We study two-dimensional Josephson arrays driven by a combined DC plus 
AC current, and with an applied transverse magnetic flux of / flux quanta per pla- 
quette. We present ansatz solutions for sufficiently large frequencies, which are a 
generalization of the travelling wave solutions found by Marino and Halsey 8 for the 
case of DC current driving. For / = 1/2 and / = 1/3, we compute the widths of 
the first few Shapiro steps for both integer and fractional winding numbers. These 
expressions consist of products of Bessel functions of (iac/^ac), where iAC and 
u>ac are the amplitude and the frequency of the driving AC current, respectively, 
times a frequency-dependent factor for fractional steps. In the limit of large fre- 
quencies, we find that the fractional steps are suppressed, whereas the maximum 
integer step widths saturate to a frequency-independent value. We show that the 
suppression of the fractional steps is due to decrease of the vertical (i.e. perpen- 
dicular to the direction of flow of the injected current) supercurrent relative to the 
normal current, whereas the persistence of the integer steps is due to the existence 
of zero- frequency (though spatially varying) terms in the expansion for the gauge- 
invariant phase differences, for which the normal current vanishes. These results 
are in reasonable agreement with the numerical computations carried out by other 
groups. 3 ' 5 
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I. INTRODUCTION 

When a resistively shunted Josephson junction is driven by a combined DC 
and AC current I(t) = Idc + I AC sin(Q Act), the current- voltage characteristic 
exhibits plateaus in which the time-averaged voltage is equal to an integer multiple 
of MlAc/^e for a finite interval of DC current. These plateaus are called Shapiro 
steps. 1 An analogous effect is observed when the current is applied to an N x N 
square array of Josephson junctions with tranverse magnetic flux per plaquette 
$ = f&a, where $o = hc/2e is the quantum flux, and / = p/q is the frustration, p 
and q being relatively prime integers. The total voltage across the array is locked 
at values given by 

nNhO, A c 

VN = -q^2e— (1) 

In this case, the steps are called fractional giant Shapiro steps. 2 The accepted 
explanation of this effect is that the q x q periodic vortex super-lattice moves 
coherently in response to the external AC field. 2 ' 3 

If the parallel shunt resistance is R, and the critical current per junction is io, 
we can then define a dimensionless time r = (2eioR/K)t. Measured in these time 
units, the external AC frequency is ujac = hVtAC /2eioR. The Josephson frequency 
is defined by fij = 2cVn /Nh, with its normalized version given by uj = Vn / iqRN . 
In terms of them, the above relation can simply be expressed as uij = vujaci where 
v = n/q is the winding number. 

These steps display a variety of characteristics. It has been observed that 
there exists a qualitative difference between the cases where v is an integer (inte- 
ger steps), and when it is a fraction (fractional steps). 2-5 This difference becomes 
manifest when going from low to high frequencies. At low frequencies, both frac- 
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tional and integer steps behave qualitatively as in the single-junction case. In the 
high-frequency limit, on the other hand, fractional steps are suppressed, whereas 
the maximum integer step widths saturate to a frequency-independent value. 4-6 

This phenomenon has been widely studied by several authors. 2-7 Analytical 
treatments have been provided by Halsey, 7 Lee and Halsey, 6 and by Rzchowski et 
al. 4 Nevertheless, a theoretical derivation based on first principles is still lacking. 
The missing link in the understanding of this problem has been the knowledge 
of the solutions for DC plus AC current driving. Progress in this direction has 
already been made. The existence of a family of travelling wave solutions for DC 
current-driven arrays has been reported by Marino and Halsey 8 in the limit of high 
Josephson frequencies. In this paper we present a generalization of these travelling 
wave states to the case where the array is driven by an additional AC current. 
For / = 1/2, these are also solutions to the model equations used by Rzchowski 
et al. 4 In addition to the modes corresponding to the Josephson frequency uj and 
its harmonics, these solutions contain terms oscillating with frequencies given by 
linear combinations of uj and ujac • If one then computes the current flowing 
across the entire array, one finds that only the terms with frequencies given by 
(kiquj + muAc) survive, where k\ and m are integers. This is due to the fact 
that these terms are exactly on phase everywhere on the array, whereas the other 
terms have phases such that their sums vanish. Shapiro steps result when this linear 
combination of frequencies is zero. We then derive expressions for the step widths 
as a function of %ac an d ojac for / = 1/2 and / = 1/3, by computing the ensuing 
DC supercurrent corresponding to this mode across the array for the first few steps. 
These expressions consist of products of q Bessel functions of (iac/mac), and other 
factors that depend solely on the frequency and an arbitrary constant phase -i/V 

The main idea ensuing from this analysis is that the difference between the 
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behaviors of the fractional and integer steps at low and high frequencies is deter- 
mined by the relative size of the vertical (perpendicular to the direction of flow of 
the injected current) supercurrent relative to the normal current for the different 
frequency modes contributing to a given step width. At low frequencies the vertical 
supercurrent is dominant for all modes, and hence, both integer and fractional steps 
behave in the same manner. When the frequency is increased, all the modes with 
non-vanishing normal currents decrease. The suppression of the fractional steps 
at high frequencies appears as a consequence of their dependence on these modes, 
whereas the persistence of the integer steps is due to the existence of zero- frequency 
modes (and therefore, with vanishing normal currents) whose amplitudes are deter- 
mined by the vertical supercurrent. Consequently, integer steps behave in the same 
fashion in both the low- and high-frequency regimes. This is the same mechanism 
that yields the Shapiro steps for a single Josephson junction. 

Our solutions for the gauge-invariant phase differences contain the arbritrary 
phase ipQ. On a given step, this phase varies over an interval that we hypothetise 
to be frequency-dependent. At high frequencies we shall assume that this interval 
attains a constant size, which in principle can only be determined by fitting the 
numerical results to the theoretical predictions. Our ignorance about the details of 
the nature of the solutions at low frequencies does not allow us to provide expressions 
for the step widths in this regime. Our study, thus, does not address the problem 
of trying to find the true variation of this interval with the frequency. 

In section II we review the solutions for the DC case and present their gener- 
alization to the case of DC plus AC current driving. In section III we use these 
solutions as starting point to compute the widths of the first few Shapiro steps for 
/ = 1/2 and / = 1/3. Finally, in section IV we provide results from numerical 
computations. 
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II. SOLUTIONS FOR DC PLUS AC CURRENT DRIVING 



We consider a square array of N x N overdamped resistively shunted Josephson 
junctions in a uniform transverse magnetic field with / flux quanta piercing each 
plaquette, parallel shunt resistance R, and critical current per junction zo- We 
define the gauge-invariant phase differences by $ij = 9i — 9j — Aij, where 9i is the 
superconducting phase on the i'th site on the array, is the line integral of the 
magnetic vector potential = ^ J! 1 A ■ dx, such that J2 P A^ = 2wf, where the 
sum is around a plaquette, and j denotes a site that is nearest neighbor with i. 
Then the current flowing from the i'th site to the j'th site is 



where 7y = Iij/io- The first term is the normal current, while the second one rep- 
resents the super current. The equations of motion simply express the fact that the 
total current arriving at each site on the array should equal the current externally 
injected there 



where the external current ii- ex t vanishes everywhere inside the array, except for 
at the boundaries. Henceforth, we shall take the convention that when computing 
the gauge-invariant phase differences, on horizontal bonds j is to be taken to the 
right of i, and above it on vertical bonds. 

For the case in which the system is driven by a uniform DC current injected 
parallel to one of the axes of the array (which we take to be the horizontal axis with 
the current flowing from right to left), and with periodic boundary conditions along 




(9ij) + sin(0y), 



(2) 




(3) 
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the vertical direction, Marino and Halsey 8 reported the existence of a family of 
travelling wave solutions. These solutions are characterized by a parameter 5 that 
measures the phase shift of the phase oscillations along the horizontal direction. 
Along the vertical direction the phase shift is simply equal to 2nf, which is consis- 
tent with the condition of transverse periodic boundary conditions with period q. 
These solutions possess a combined spatio-temporal translational symmetry, in the 
sense that a translation of the solution by one lattice spacing along the horizontal 
direction is equivalent to the translation of the solution by a time r = S/uj, and 
a translation of the solution by one lattice spacing along the vertical direction is 
equivalent to a translation of the solution by a time r = 2irf/ujj. 
Let us define 

ip = ujt + 2tv fn Y + Sn x + i>o, (4) 

where nx and ny are integers, and tpo is a constant phase. Then, the gauge- 
invariant phase differences on horizontal and vertical bonds for these solutions are 
given by 

MV0 = V> + MV0 ; (5) 
MVO = MV0> (6) 

where fu and fy are periodic functions with period 2tt and zero time average. The 
authors of reference 8 worked out the analytical form of these functions in the limit 
of high voltages, by retaining only the first harmonic in their Fourier expansion. 

The generalization of these solutions to the case of combined DC and AC 
current driving is straightforward. The main modification is that now in addition 
to the mode with Josephson frequency uj, a mode with frequency ujac should also 
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be present, due to the external AC current. The presence of this second time scale 
ruins the spatio-temporal translational symmetry in the cases in which ui j and ujac 
are incommensurate. This entails no problem, as we shall see. The beating of these 
two modes due the horizontal supercurrent requires the additional presence of terms 
with frequencies given by linear combinations of uij and ujac in the expansion for 
the gauge-invariant phase differences. These solutions must further satisfy several 
conditions. Firstly, they should reduce to their DC counterparts when %ac and 
ujac are set equal to zero; secondly, in order to (indirectly) enforce the boundary 
conditions, the DC and AC currents (at ujac and its harmonics) flowing on each 
horizontal bond should be independent of position, and equal to the values given 
by the currents injected at the boundaries, whereas on vertical bonds they should 
vanish, allowing for the possible existence of zero- frequency modes (that occur only 
on integer Shapiro steps), that should not be regarded strictly as DC currents; and 
finally, the linear term should remain unchanged, because it is still true that the 
slope should yield the average voltage per junction. Our ansatz for the gauge- 
invariant phase differences on horizontal bonds then takes the form 

H (n x , n Y ,r) =i) + Y^ ®n, m cos(mf) + muj A c r + f£ m ), (7) 

n,m 

where n (> 0) and m are integers, while on vertical bonds we have 

e v (n x ,n Y ,r) = ^ o^ m cos(ra/> + moj AC r + f£ m ), (8) 

n=£0,m 

where ip is given by Eq. (4), and are constant phases. The phase differences 
are taken according to our convention. The different terms in this expansion are 
labeled by the two integers n and m. We shall refer to this component of the 
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phases and currents as the (n, m) mode. Notice that the modes with n = have 
no spatial dependence, in agreement with our assumption. This form of solution is 
good enough to describe the system even at low frequencies. 

The gauge-invariant phase differences are not independent. The sum of their 
oscillating parts around a plaquette has to vanish. We impose this condition to 
each frequency mode and obtain 



where (3 = sinnirf / sin(n5/2), if n ^ q (n = q is a shorthand for n = q mod 0), 
and otherwise. Thus, the modes with n = q are absent on vertical bonds. This 
result is the same one that was obtained in the DC case. The reader is refered to 
reference 8 for the details of the derivation. 

On the Shapiro steps we will assume that 6 = 2nf, due to our requirement of 
AC translational invariance (see the discussion preceeding Eq. (7)). This implies 
(3 = 1 for n 7^ q, in which case a and £ are the same on both horizontal and vertical 
bonds, which considerably simplifies matters. Consequently, we shall hereafter drop 
the superscripts H and V in our expressions. 

We shall perform a mode expansion of the horizontal supercurrent sm9n in 
the following manner: 



a. 



y 

n,m 




(9) 



Cm = Cn,m + n(7Tf-S/2), 



(10) 




n,m 



(11) 




n,m 
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The different components of the supercurrent can be computed in a straight- 
forward manner by performing a Fourier-Bessel expansion of the left-hand side of 
Eq. (11). A generic term in this expansion is of the form 




j 

The sum is over the set of sets of pairs {{(nj,mj)}} such that there exists a 
set of coefficients {kj} for which the following relationship (understood as a vector 
identity) holds 

(n,m) = (l,0) L + ^2 kj(nj,mj), (13) 
j 

where (1,0)^ denotes the contribution due to the linear term, which is absent on 
vertical bonds. These expressions are in general quite complicated. We shall assume 
that they can approximately be computed starting from the lowest-order (in n and 
m) modes, since the magnitude of the different Bessel functions decay quickly with 
increasing order. This approximation should be good enough at high frequencies, 
but we do not expect it to remain accurate for lower frequencies. The feedback of 
the higher-order modes on the lower ones should become more important as one 
approaches the critical current. On vertical bonds the (n, m) component of the 
supercurrent is (neglecting higher-order corrections) equal to 2Ji(a n m ). It can be 
shown for simple cases that other contributions vanish. 

We now turn to the equations of current conservation (Eq. (3)). Once again, 
we have to distinguish between the cases n ^ q and n = q . In the former case, the 
equation for current conservation for our ansatz solution is 
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-2u ntrn a n , m sm(nip + mu AC r + f n>m ) + 2Ji(a n>m ) cos(mp + mu AC r + £ n , m ) 

+ S njm cos(ntp + mu AC T + E njm ) = 0, (14) 

where u; n)m = nwj + mu A c- There is an overall factor of 2siri7r/ that goes away. 
The first term in the above equation represents the combined effect of both the 
horizontal and vertical normal currents, each of them contributing the same amount; 
the second one is due to the vertical supercurrent, and the last one comes from the 
horizontal supercurrent, which has to be computed for each mode. For n = q the 
current is trivially conserved at each site. Furthermore, the current corresponding 
to the (0, 1) mode should equal the external AC current: 

-u AC a ,i sm(u AC t + f 0) i) + S 0A cos(u AC r + S ,i) = -i A c sin(u AC r). (15) 

For convenience, and without loss of generality, we have introduced a minus sign 
at the right hand side of this equation. If we neglect the vertical supercurrent and 
replace the factor of 2 multiplying the first term in Eq. (14) by 1, then these two 
last equations describe an overdamped single junction. 

In general, we expect cto,i ~ ^ac/^ac^ an d the asymptotic behavior of the 
different amplitudes at high and low frequencies to be of the form 

h{uj) n j£(«o,i)- (i6) 

rij ki=m 

for n 7^ q, where h(x) ~ c (c = constant < 1), as x — > 0, and h(x) ~ x~ n , as 
x — > oo. This will be made clear below. Consequently, it is safe to assume (save 
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for ao,i) that Ji(a n ,m) ~ «n,m/2. Eq. (14) can then be solved in terms of the 
expression for the horizontal supercurrent: 

Sn,m (17\ 

OLn,m — i i \*- I ) 

^nm = 7T + S n , m - arctan(2w n>m ). (18) 

For uJn^m 3> 1 the effect of the vertical supercurrent can be neglected and 
«n,m ~ 5 , n,m/(2cUn,m) 5 whereas for u n ^ m <C 1, the vertical supercurrent dominates 
and a n ^m ~ S niTn . Eq. (16) can then be proven in the following manner. Since 
(0, m) = (l,0)i + m (0,1), then, to leading order So }7n = Jm^Ac/^Ac)- From 
Eq. (17) it follows that «i jm has the asymptotic behavior given by Eq. (16). The 
proof for OL n ^m can be made by induction. 

In the present calculation we will neglect the supercurrent in Eq. (15). This 
is a good approximation for ujac 1. Thence, cto,i = ^ac/^ac and £o,i = 0. 
Using this, (1, 0) = (1,0)^ + (0, 1), and Jo(an,m) ~ 1 for all the other modes, we 
find 



n 



Ci.o = g ~~ arctan(2wj). (20) 

This solution also holds in the DC case (zac = 0), and represents a generalization 
of the solutions presented in reference 8. Unlike them, this solution remains regular 
as ujj — > 0, owing to the vertical supercurrent. Similarly, (1, ±1) = (1, 0)l ± (0, 1). 
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Then, 



«i,±i = 



Jl(iAc/vAc) 



(21) 



^4(UJJ±LVAC) 2 + l' 



£i,±i = 7r - arctan(2(wj ± u AC )). 



(22) 



Other modes can be computed in a similar fashion. 

III. THE SHAPIRO STEPS 

It is clear that both the normal current and the supercurrent have the same 
harmonic dependence as the gauge-invariant phase differences. In particular, this 
implies that when computing the total voltage and supercurrent across the array, 
only the modes with n = k\q survive. It is immediate to check that all the other 
terms cancel out. This is usually interpreted in terms of the vortex configuration 
by saying that a vortex moves q times during a period of 2tt/ujj. This is key in 
order to understand the phenomenon of the Shapiro steps. 

Looking back at Eq. (11), in presence of an external AC current, an additional 
DC supercurrent will appear across the array for frequencies satisfying k\qujj + 
muJAc = 0, with ki and m relative primes. The Shapiro steps ensue. We see 
that the proposed rigid motion of the q x q vortex super lattice in response to the 
external AC field has a very natural explanation within our theory. For k\ > 1 we 
have subharmonic steps. We will not consider this possibility here, because these 
steps are in general too small to be observed. 

The first step in our calculation is to identify the modes that yield the largest 
contribution to the DC supercurrent corresponding to the different Shapiro steps. 
For integer steps {v = n) the choice is unambiguos: it is the set of zero-frequency 
modes (&;', — k'n), where 1 < k! < q. These modes have the virtue that their 
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associated normal currents vanish, and thus, they are determined by the vertical 
supercurrent. The Shapiro step widths in this case depend only on (iac/^ac)- The 
case of fractional steps {y = n/q) can be analyzed in a similar manner. The most 
important contributions are due to the modes (n',n"), with n' < n and n" < q. 
The amplitudes of these modes decay with the frequency because they are mainly 
determined by the normal current. 

The distinction between the low- and high-frequency behaviors just amounts to 
saying that at low frequencies both fractional and integer steps are in a supercurrent- 
dominated regime, whereas at high frequencies only the integer steps are, thanks to 
the zero-frequency modes. Saying that the steps display single-junction behavior is 
just another way of rephrasing this fact. 

All of our expressions depend on the arbitrary phase ipo. In particular, the 
DC supercurrent on a given step depends on this phase, and thus the width of 
the step will depend on the range of variation of it. The interval of variation of 
this phase should depend on the dynamical stability of these solutions and also 
on the nature of the solution for the (0, 1) mode at low frequencies (recall that 
we neglected the supercurrent term in Eq. (15)). We conjecture that the size of 
this interval varies with the frequency, growing from zero to an interval of constant 
size at high frequencies, and that this is the mechanism underlying the growth and 
saturation of the steps. This assumption seems to be good enough to reproduce 
all the observed qualitative features of the steps. The quantitative correctedness of 
our expressions will depend on the goodness of our guess for the variation of this 
phase. According to this, our results for a given frequency should differ at most by a 
constant factor (for all values of the iac) from the results obtained from simulations 
or experiments. We shall not attempt to resolve this issue here, rather, we shall use 
a generalization of an ansatz used by Halsey, 7 that seems to work well to reproduce 
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the observed numerical results in a reasonable manner. Our results have a factor 
of the form cos (qipo + 4>q,u)- We shall assume that at high frequencies (qifro + 4>q,u) 
is centered at tv/2 for q even, and at for q odd. Furthermore, for both q even 
and odd, we shall assume that this phase varies over the interval [— 7r/2,7r/2] for 
fractional steps, and over [— n/2q 7 ir/2q] for integer steps. We shall thus restrict 
ourselves to making predictions for the values of the step widths only for large 
enough frequencies. 

We now turn to specific examples. We shall only consider the cases / = 1/2 
and / = 1/3, because the number of modes that have to be included in a given 
calculation increases quickly with q. 

a. Integer Steps: v = n 

(al) / = 1/2 

A DC supercurrent occurs at the mode (2, — 2n) = (1, 0)l + (1, —n) — n (0, 1). 
In this case «i,- n = Jn(iAc/<^Ac), and £i,- n = (n + l)ir/2. We thus obtain the 
expression: 

k/2,n = Ji(ai,-n) 4(«o,i) sin(2^ + (2n + l)vr/2) 

= g {Jn{iAc/ooAc)f sin(2V> + (2n + 1)tt/2), (23) 

where i\/2,n denotes the ensuing DC supercurrent. This is not yet the expression 
for the step width. The final answer depends on the range of variation of ip . Using 
our ansatz for this variation at high frequencies we find 



Azi/ 2 ,n = {Jn(}Ac/uAc)) 2 - 



(24) 
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(a2) / = 1/3 

Here we have two contributions, namely, (3, — 3n) = (1, 0) l + (2, —2n)—n (0, 1), 
and (3,-3n) = (1, 0) L + 2 (1, -n) - n (0, 1). Now, a 2 ,_ 2n = ( J n {i A cl^Ac)? 1^ 
and ^2,-2n = (n + 3/2) 7r. The total DC supercurrent corresponding to this step is 

H/3,n = Jn(iAc/uAc) -Mai -n) sin(3^ + 3mr/2) 

+ Ji{a-2,-2n)Jn{iAc/uAc) sin(3-i/>o + 3n7r/2) 
= \ {Jn{iAc/oo AC )f sin(3^ + 3mr/2). (25) 

o 

At high frequencies this becomes 

Aii/ 3l „ = ^(4(ucMc)) 3 (26) 

In conclusion, we find that the integer step widths are independent of the frequency. 
It should be clear that in general i p / q ,n ~ (Jn(iAc/ UJ Ac)) q ■ F° r 9=1) this reduces 
to the single-junction result, or equivalently, to the result for the unfrustrated case 

f = o. 

b. Fractional Steps: v = n/q 

We need to compute the ((/, —n) component of the horizontal supercurrent. As 
mentioned above, the most important contributions come from the modes (n', n") 
such that n' < q and n" < n. The number of modes that are relevant for the 
determination of the steps grows quickly with n, though, so we will only work out 
explicitly a few simple cases. 
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(bl.l) / = 1/2, 1/= 1/2 

We can either write (2, -1) = (1, Q) L + (1, -1) + (0, 1), or (2, -1) = (1, 0) L + 
(1, 0) — (0, 1). Omitting some details we get 



H/2,i/2 = - Jo(oto,i)Ji(oti,-i) cos(2V> + arctan^Ac) 
- ^i(ai,o)«/i(Q!o,i) cos(2^ - arctancuAc) 
Using the expressions for 0:1,0, ^o,!, an d 01,-1, this turns into 



(27) 



H/2,1/2 = 7-2 ~ ^ cos(2^o). (28) 

and at high frequencies 

a- 2Jo(iAc/uAc)Ji(iAc/uAc) /on x 

Az V2,l/2 = 7-2 7^ ( 29 ) 

We find that the step width decays like l/u) 2 AC at high frequencies, in dis- 
agreement with what has been assumed by other authors. 4 ' 5 At low frequencies, 
this result reduces to a frequency-independent expression, which is characteristic 
of single-junction behavior. It is also possible to compute higher-order correc- 
tions (in Bessel functions of (iac/^ac)) by considering the combinations of modes 
(2,-1) = (1,0) L + (1,1) - 2(0,1) = (1,0) L + (1,-2) + (0,1). It is not hard to 
see that these corrections go like Ji{Iac /^ac) J2(iAc/ UJ Ac)/(^ UJ AC + -0- This is 
negligible compared to the expression given in Eq. (29) for most cases of interest. 

(bl.2) / = 1/2, v = 3/2 

In this case we have (2, -3) = (1, 0) L + (1, -l)-2 (0, 1), and (2, -3) = (1, 0) L + 
(1, —2) + (1, —1) + (0, 1). We will neglect other contributions. The calculation is 
identical to the previous one and yields 
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H/2,3/2 = 7-2 ,^ COS(2^ )- (30) 

In the limit of high frequencies 

1/2,3/2 = (^TTT) (31) 

The next-order correction varies like Jq^ac/^ac) J^ac/mac), which can be ne- 
glected. We see that this step width has the same dependence on ujac as for 
v = 1/2. 

(b2.1) / = 1/3, v = 1/3 

This calculation involves higher-order modes than the ones hitherto used . The 
number of terms, therefore, considerably increases. In fact, we now have the possible 
combinations: (3, -1) = (1, 0) L + (1, -1) + (1, 0) + (0, 1), (3,-1) = (1,0) L + 
2(1,0)-(0,1), (3,-1) = (l,0) L + (2,0)-(0,l),and (3, -1) = (1, 0) L + (2, -1) + 
(0, 1). There are two contributions to the (2, —1) mode, which will be considered 
separately. These are (2,-1) = (1, 0) L + (1, -1) + (0, 1) = (1, 0) L + (1, 0) - (0, 1). 
Putting everything together we obtain: 



Jo(iAC'/uAc)Ji(iAc/uAc) ( cos(3V>o - 2 arctan(2^Ac/3)) 
H/3,i/3 = — 



+ 



V4u^ c /9 + l V 2^4^/9 + 1 

cos(3V'o + arctan(2o; y ic/3) + arctan(4u; J 4c , /3)) 
a/16^/9 + 1 



cos(3V'o + arctan(4o;Ac/3) — arctan(2o;/ic , /3)) 
+ v/16^/9 + 1 (32) 

cos(3-i/>o — arctan(4o; y ic/3) — arctan(2u;,4c73)) 



+ 



v/16^/9 + 1 

cos(3V'o) 



+ 1 
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(33) 



After some algebraic manipulations this becomes 



_ Jl{i A c/u A c)Ji(iAc/uAc) \( 3-8c4 c /9 , 2^ c /9 + 3/2 \ 
H/3 ' 1/3 4(4^ c /9+l) "[^(16^ c /9 + l) + (4^/9 + 1)^ 

/<j / \ i ( 2uac/3 2uj ac /3 \ . . ' 

COS(3 "'°) + ( (M c /9 + l) - (16^/9 + d J Sm(3 * ) 

In the limit of high frequencies this is 
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^1/3,1/3 = 10g 3 j o(*acMc) J\{iAc/wAc) sin(3^o), (34) 
and with our ansatz for 

81 

Aii/3,1/3 = 1oa 3 Jo(iAc/wAc) Ji(}ac/wac) (35) 
whereas in the low-frequency limit we find 
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H/3,1/3 = g Jo{iAc/u A c)Jl{iAc/u A c) COs(3V>o)- (36) 

(b2.2) / = l/3,i/ = 2/3 

The number of modes to be included in the calculation keeps growing, as 
promised: (3,-2) = (1, 0) L + 2 (1, -1) + (0, 1), (3,-2) = (1, 0) L + 2 (1, 0) - 
2 (0, 1), (3, -2) = (1, 0) L + (2, -1) - (0, 1), (3, -2) = (1, 0) L + (2, 0) - 2 (0, 1), 
(3, -2) = (1, 0) L + (2, -2) + (0, 1), (3, -2) = (1, 0) L + (1, -1) + (1, 0) - (0, 1), 
and (3, -2) = (1, Q) L + (1, -2) + (1, 0) + (0, 1). In order to carry out the calcula- 
tion we need (1,-2) = (1,0)^-2(0,1), and (2, -2) = (1, 0) L + (1, -1) - (0, 1) = 
(1, 0) L + (1, 0) - 2 (0, 1) = (1, 0) L + (1, -2) + (0, 1) . The calculation is analogous 
to the one done in the previous section. In the high-frequency limit we obtain: 
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81 



Jo(iAcfuAc) {Ji(iAc/uAc) 
+ Jo(iAc/u A c)J2(iAc/uAc)/8) COs(3^o), 



Zl / 3 ' 2 / 3 " 128^! 



AC 



(37) 



and after the usual assumption for ipo 



Ai 



81 



Jo(iAc/wAc) {Ji(iAc/uAc) 
+ Jo(iAc/uAc)J2(iAc/u A c)/&)- 



1/3,2/3 



128 u\ 



AC 



(38) 



We have included higher-order corrections in Bessel functions that this time yield 
a no n- negligible contribution. Once again, the step width decays as l/u AC . The 
low-frequency limit yields 



We get the expected single-junction behavior. 

We see that at high frequencies the step- widths for / = 1/3 decrease like l/u> AC . 
The terms that vary like 1/uj 2 ac cancel out. Compare this result with the ones ob- 



The term varying like 1/uac also did cancel out. This result seems to be quite 
general, and there is a way of understanding it. The fractional steps in the present 
case correspond to the subharmonic steps for a single junction. We can say that 
the vertical supercurrent plays the role of an effective additional degree of freedom, 
analogous to the role played by the capacitive term for a single junction, in which 
case subharmonic steps do appear. At high frequencies the vertical supercurrent 



H/3,2/3 = 



9Jp(iAc/uAc) 
8 




tained for / = 1/2. In those cases we found that the step width varied like l/u AC . 
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becomes negligible, and then the equation for current conservation effectively be- 
comes the equation for an overdamped single junction, for which the subharmonic 
steps are absent. 

Another observation that can be made at this point is that in all the cases that 
have been studied, the expressions of the step widths for / = p/q and v = n/q 
involve the products of q Bessel functions of (iac/^ac) such that the sum or the 
difference of their orders is equal to n. This is in accord with Eq. (16). The pre- 
factor has a different dependence on the frequency due to cancellations occuring 
among the different modes. 

IV. NUMERICAL RESULTS 

We use the same method of numerical integration used in reference 8. We 
briefly review it here, for convenience. The equations for current conservation can 
be written as a matrix equation 

rift 

My^ = F({0,-<W}) (40) 

where %' denotes a nearest neighbor to i. The matrix M is then inverted yield- 
ing a set of coupled first-order differential equations which we integrate using the 
fourth-order Runge-Kutta method. 9 The current was uniformly injected at the left 
boundary of an N x N array. Furthermore, we used periodic (with period q) bound- 
ary conditions in the direction perpendicular to that of injected the current. We 
normally used arrays of size N = 3q to 5q, in order to avoid effects due to the 
boundary conditions. We used staircase configurations as initial states in all of our 
simulations. We restricted our observations to the cases / = 1/2 and / = 1/3, for 
the already mentioned reasons. 
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Fig. 1 shows the power spectrum of the oscillating pieces corresponding to the 
gauge-invariant phase differences on horizontal and vertical bonds for / = 1/2, away 
from any step. The slope of the linear term on horizontal bonds agrees with the value 
of the Josephson frequency observed in the power spectrum. Peaks 3 though 6 are 
absent on vertical bonds, in particular notice the absence of the (0, 1) mode. This 
is accord with our assumptions. This same behavior has been observed for / = 1/3. 
Peak 4 in Fig. 1 (a) corresponds to the (2, —1) mode. On the first fractional step the 
frequency corresponding to this mode is zero, and the corresponding supercurrent 
yields the additional DC current on the step; the same is true for Peak 7 on the first 
integer step. Also, on the latter step, the mode corresponding to Peak 1 corresponds 
to the zero frequency mode. On the different steps we get the same picture for the 
power spectrum, with the difference that in these cases the motion is periodic. 

Fig. 2 displays the behavior of the gauge-invariant phase differences on hori- 
zontal and vertical bonds for / = 1/2 and v = 1. The presence of a zero- frequency 
(spatially varying DC component of the phase) is clear. This is also in agreement 
with our results. 

In Fig. 3 we show the different step widths as a function of %ac for ujac = 2, 
and ujac = 3. At higher frequencies the fractional steps become suppressed. These 
results are in agreement with the simulational results obtained by other groups. 3 ' 5 

In Fig. 4 we show the variation of the step widths with %ac for ujac = 2.0, 
and ujac = 3.0. The qualitative behavior of the step widths is the same that was 
found for / = 1/2. 
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FIGURE CAPTIONS 

1. (a) Power spectrum P(u) of horizontal phase oscillations for / = 1/2 away from 
any step. The peaks (1 — 6) correspond to the frequencies (ujac — ujj), ujj, 
ujac, (2cuj — ujac), 2ujj, and (Iujac — 2cuj) respectively; (b) Power spectrum of 
vertical phase oscillations. Notice the absence of peaks 3-6, which is in accord 
with the selection rule for the existence of modes on vertical bonds (see the 
remark following Eq. (10)). 

2. Gauge-invariant phase differences on two vertical bonds lying on the same 
column for / = 1/2 on the first integer step. They are out of phase by S = n, 
as assumed in our work. The average value of the phase oscillations is non- 
vanishing and is depicted by a solid line in the figures. This represents the 
zero- frequency mode, which changes phase by tt when translated by one lattice 
spacing on the array. 

3. Step widths for / = 1/2 and v = 1/2, 1, 3/2, and 2 for (a) ujac = 2.0; (b) and 
ujac = 3.0. Fractional steps are suppressed relative to the integer steps. 

4. (a) Step widths for / = 1/3, v = 1/3, 2/3, 1 , and ujac = 2.0; (b) The 
frequency is ujac = 3.0. The steps are considerably smaller than those for 
/ = 1/2- 



